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We numerically study imbalanced two component Fermi gases with attractive interactions in 
highly elongated harmonic traps. An accurate parametrization formula for the ground state energy 
is presented for a spin-polarized attractive Gaudin-Yang model. Our studies are based on an accurate 
microscopic spin-density-functional theory through the Kohn-Sham scheme which employs the one- 
dimensional homogeneous Gaudin-Yang model with Luther-Emery-liquid ground-state correlation 
as a reference system. A Thomas-Fermi approximation is examined incorporating the exchange- 
correlation interaction. By studying the charge and spin density profiles of the system based on these 
methods, we gain a quantitative understanding of the role of attractive interactions and polarization 
on the formation of the two-shell structure, with the coexisted Fulde-Ferrell-Larkin-Ovchinnikov- 
type phase in the center of the trap and either the BCS superfluid phase or the normal phase at the 
edges of the trap. Our results are in good agreement with the recent theoretical consequences. 

PACS numbers: 03.75.Ss,71.15.Mb,71.10.Pm 



I. INTRODUCTION 

Progress made in trapping ultracold atomic gases into 
different low-dimensional systems [l|, Q has provided us 
grounds to describe and simulate the behavior of fasci- 
nating non-Fermi liquid in which the analysis may not 
be as tractable. The physical properties of such systems 
have attracted both experimental and theoretical inter- 
est. 

For the research in the experimental side, the Tonks- 
Girardeau regime of strongly repulsive bosons in one- 
dimensional (ID) systems has been observed [2]. This 
regime is understood through the mapping of the strongly 
repulsive, impenetrable bosons onto an ideal gas of 
fermions subjected to the same external potential 
Esslinger and colleagues at ETH Zurich in Switzerland 
confined a two-component Fermi gas of 40 K atoms to 
thousands of highly elongated one-dimensional tubes. 
Such a system is realized by using a two-dimensional op- 
tical lattice, having about 100 atoms in each tube |l| 
and serves as a ID matter waveguide to observe the 
two-particle bound states of atoms. For the first time, 
they created confinement induced molecules. The bind- 
ing energy of the molecules is measured by using radio- 
frequency spectroscopy bearing a good agreement with 
the theoretical predictions Q. Importantly, two-particle 
bound molecular states are formed in ID confined system 
independent of the nature of the interaction, i.e., attrac- 



tive or repulsive between atoms. This is obviously differ- 
ent from what happens in free space where a molecule is 
formed only when the atoms attract. 

On the other hand, many theoretical efforts on quasi- 
one dimensional (Q1D) inhomogeneous Fermi gases have 
been made. For ID trapped fermions at finite tem- 
perature with attractive binary interactions, the tran- 
sition to a fermion-paired state is obtained by exact 
stochastic mean- field calculations 0, which is charac- 
terized by dominant Cooper pairing correlations, an al- 
gebraic long-range order and a superfluid component. A 
one-dimensional spin-polarized Fermi gas with infinitely 
strong attractive zero-range odd-wave interactions has 
been studied. This system is known as the fermionic 
Tonks-Girardeau gas, arising from a confinement-induced 
resonance reachable via a three-dimensional p-wave Fes- 
hbach resonance In Ref. 0], a bosonization tech- 
nique has been developed to calculate analytically the 
density profile, the momentum distribution, and sev- 
eral correlation functions of two-component Fermi gases 
with inclusion of forward-scattering processes and exclu- 
sion of backward scattering between inter-components. 
Inhomogeneous Tomonaga-Luttinger liquid model, with 
space-dependent parameters assuming a slowly varying 
trap potential on the scale of the Fermi wavelength, 
has been used to describe spin-charge separation in two- 
component Fermi gases 0, Q. The Thomas- Fermi ap- 
proximation is widely used to calculate the ground state 
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properties of a large system 0, [l^|. The existence of 
a molecular Tonks- Girardeau gas is predicted for strong 
attractive effective interactions [10]. In Ref. 11 1 a mi- 
croscopic calculation of the density-functional theory 
(DFT), based on the exact Bethe-ansatz solution of ID 
homogeneous Gaudin-Yang system, is presented for the 
ground-state density profile with arbitrary size. Quan- 
titative understanding about the role of the repulsive or 
attractive interactions on the shell structure of the axial 
density profile is achieved. 

Above mentioned studies are focused on the fully po- 
larized Fermi atomic gases. Recently, progress in two ex- 
periments [3, E3] in trapping partially polarized Fermi 
atomic gases has paved the way for systematic studies on 
systems with different number of spin-up and spin-down 
atoms (ATf ^ iVjJ. The experiment explored on three- 
dimensional polarized 6 Li gases 12j suggested that the 
density profiles show coexisting three-shell structures: a 
fully unpolarized paired superfluid phase in the center 
of the trap, a fully polarized noninteracting phase com- 
posed of excess atoms, and a partially polarized normal 
shell between these two regions. The experiment which 
employed an elongated cigar-shaped trap of polarized 6 Li 
gases [l3[ suggests the gas separates into a phase that is 
consistent with a superfluid paired core surrounded by a 
shell of normal unpaired fermions beyond a critical po- 
larization. 

Recently two theoretical investigations El dis- 
cussed partially polarized two-component Fermi gases of 
attractive interaction within local density approximation 
while implementing the exact exchange-correlation inter- 
action from the homogeneous Bethe-ansatz solution of 
ID Gaudin-Yang model. Both of these two papers predict 
a two-shell density profile structure of a coexisting par- 
tially polarized superfluid core in the center of the trap 
while at the edges of the trap either fully paired or fully 
polarized wings depending on the attractive interaction 
strength and the polarization. Furthermore, they antic- 
ipated a critical spin polarization and a non-monotonic 
behavior of the Thomas-Fermi radius of each spin com- 
ponent as a function of the polarization. Orso L| and 
Hu et al. [l5| identified the polarized superfluid as hav- 
ing the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) struc- 
ture, which is stabilized in a polarized two-component 
Fermi gas in an array of weakly-coupled ID tubes 16 1. 
These two theoretical efforts used the local-density ap- 
proximation for the harmonic trap and Thomas-Fermi- 



type approximation for the interaction. The local-density 
approximation is applicable when TV 3> 1, where N is 
the number of particles, but deteriorates with decreas- 
ing particle number and completely neglects the tun- 
neling of the density profile beyond the Thomas-Fermi 
radius where one instead needs to use a method be- 
yond the local-density approximation, like the micro- 
scopic density-functional theory. 

The FFLO phenomena of imbalanced fermion popu- 
lations in ID optical lattices of attractive interactions 
have also been extensively discussed recently 0, 0, Q , 
where FFLO pairing is formed robust from the attractive 
interaction for a wide range of polarization, and for the 
spatial inhomogeneities induced by the presence of a trap 
in the experimental reality. 

In this paper, we study ID harmonically trapped Fermi 
gases of imbalanced spin population with two alternative 
approaches. One is the spin-density- functional theory 
based on the Bethe-ansatz results for the ID Luttinger 
liquid, which is specially useful for systems whose an- 
alytical solutions or exact numerical methods are not 
available, and the other is the Thomas-Fermi approxi- 
mation based on the non-interacting approximation for 
the kinetic energy. Both methods consider exactly the 
exchange-correlation energy of the corresponding ID ho- 
mogeneous system. Within these two methods, we quan- 
titatively study the systems of weak and intermediate 
coupling strength and different polarizations by analyz- 
ing the charge and spin density profiles. 

The paper is organized as follows. In Sec. II we intro- 
duce our model and present the parametrization result 
for the ground state energy of the homogeneous Gaudin- 
Yang model. In Sec. Ill the spin-density-functional the- 
ory is briefly discussed. Numerical results and discussion 
are presented in Sec. IV. Finally, we conclude in Sec. V. 



II. THE MODEL AND ITS 
PARAMETRIZATION RESULT 



A system of N spin-1/2 fermions is considered to be 
one-dimensional when subjected to a strongly anisotropic 
harmonic potential with angular frequencies in the radial 
direction u±_ much larger than that in the axial direction 
a; 1 1 with u>n/u)±_ <C iV _1 . We consider two kind of fermion 
particles with mass m, interacting via a <5-type contact 
with an effective ID coupling strength in the har- 
monic trap. The Hamiltonian for such a system is given 
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by 



H 



h 2 N d 2 N1 Nl 
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Note that when p-wave interactions are neglected, the 
contact s-wave interactions between same species are sup- 
pressed by the Pauli exclusion principle. The coupling 
constant g\r> can be written in terms of the scattering 
strength as gio = —2h 2 /(aiDm). The effective fD scat- 
tering length a\r> can be expressed through the three- 
dimensional scattering length a^D for fermions confined 
in a quasi- ID geometry, clid — — aj_(l — Aa%D / / > 
with the constant A ps f.0326 0]. In the follow- 
ing, we will choose the harmonic-oscillator length an = 
yjh/(muj\\) as unit of length and the harmonic-oscillator 
quantum fkJn as unit of energy. 

Without the external potential the Hamiltonian (fTJ) re- 
duces to the homogeneous Gaudin-Yang model, which is 
exactly solvable by use of the Bethe-ansatz method [2o| . 
In the thermodynamic limit, the system is determined 
by the spin polarization £ = (iV-j — Ni)/N and a dimen- 
sionless parameter 7 = mgi^,/(h 2 n), while in the fully 
unpolarized system the only parameter of the system is 
7. We denote the charge density as n = N/L and the 
spin density s = (A7 - N{)/2L = (n T - nj.)/2. 

For N attractive fermions, the ground state is de- 
scribed by the coupled integral equations for the momen- 
tum distribution p{k) and <r(A), 



p(k) 
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and 
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P (k) 



(>yn) 2 + (X-Xi) 2a{X ' h (3) 



where B and Q are non-negative numbers related to the 
normalization conditions 



and 



dkp(k) = 2s 



dXa(X) = — — s . 
> _ 
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FIG. 1: (color online) The interaction contribution f((,x) = 
[£gs(?i, Ci 7) — n(n, Q]/Ef to the ground state energy of the 
homogeneous Gaudin-Yang model as a function of the spin 
polarization £ for (a) 7 = —2 and (b) 7 = —20. The ex- 
act result, obtained from the solution of the Bethe-Ansatz 
Eqs. ©-©I is compared with the parametrization formula 
given by Eq. (0. 



The ground state energy (GSE) per particle is written in 
terms of the momentum distributions p{k) and tr(A) as 



E 
N 



n 2m 



+2 



dkk 2 p(k) 



dX [A 2 - ( 7 n/2)V(A) 



(6) 



For N repulsive fermions interacting via a <5-type con- 
tact, the correlation energy per particle, and then the 
GSE, has been parameterized based on Pade approx- 

0, Si] 



imant in 



21 



Magyar achieved an improved 
parametrization results for the correlation energy by in- 
cluding the higher-order correlation kernel terms at the 
high-density limits which is needed for the adiabatic 
time-dependent DFT approximation 23 1. 

For the present system of attractive contact interact- 
ing, we solve numerically a set of coupled Eqs.j2])-([5]) 
and find the GSE per particle £gs(^, C>7) = E/N given 
in Eq. |B| as a function of n, £ and 7. After having 



4 



£Gs( n ' t y)f E F 




FIG. 2: (color online) The ground state energy per particle 
in unit of Fermi energy as functions of £ and 7. 



£Gs(n, Cj 7)) we find an accurate parametrization formula 
to fit £gs(^i Cj 7)- Our parametrization formula is 



ecs(n, C, 7) = k(u, () + /(C, x)£f 



where 



K(n,C) = 



24m 



-(1 + 3C 2 ) 



(7) 



(8) 



is the kinetic energy of the noninteracting system per par- 
ticle and the Fermi energy is E-p = h 2 n 2 ir 2 / (8m). Here, 
/(£, x) can be represented by, 

/(Cs) = [e(x)-l/3]{l + a(x)|C|+/3(x)C 2 

- [l + a(x)+/?(x)]|C| 3 }, (9) 



where x = 27/71-, e{x) is given in Ref. very accurately 
as, 

,s _ 1 \x\ x 2 + a m \x\ + b m x 2 
3 7T a; 2 + Cm|ar| + d m 4 

with a m = -0.331117, 6 m = 0.458183, c m = a m + 4/vr, 
andd m = 4a m /TT+b m +16/n 2 — 1. Equation (fT0|) gives the 
exact asymptotic behaviors at x — ■> _ and x — > — 00 [lF 
a(x) and /3(x) are given by Pade approximant, 



a(x) 
0(x) 



A a x 



- B a x + C a 
Apx + Bp 



(11) 



x 3 + Cpx 2 + Dpx - Bp 



Here A a = -0.0652385, B a = -1.87353, C a = 3.08873, 
A fJ = 3.90515, Bp = 19.4239, Cp = -4.53457 and Dp = 
-7.29826. 
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FIG. 3: (color online) The magnetization (in units of the in- 
verse of arbitrary length a) as a function of magnetic field (in 
units of h 2 /ma 2 ) for different fillings at n = 0.5 and n = 1.0 
for given 7 = — 10 (n is taken also in units of the inverse 
of arbitrary length a, in the text we take the same units as 
described here.). The exact results from Bethe-Ansatz equa- 
tions are compared with those derived from the parametriza- 
tion formula. 



To assess the validity of our parameterized formula, we 
show in Fig. [1] the results emerging from Eq. |(7J) in com- 
parison with the exact Bethe-Ansatz results for 7 = — 2 
and 7 = —20 in the weak and strong attractive regime, 
respectively. As it is clear from Fig. [T] our parametriza- 
tion formula has an excellent agreement with the exact 
calculation and it would be safe to use the parametriza- 
tion formula for all range of n, £ and 7 values. 

Here we would like to deal with some asymptotic be- 
haviors of the simple parametrization formula in Eq. (f7|). 
First of all, we can simply obtain both unpolarized 
(C = 0) and polarized (£ = 1) limits exactly. Secondly, 
the weak interacting regime [24J is recovered after taking 
7^0 limit in Eq. as 



£Gs(n,C,7 -> 0) 



fi 2 ?! 2 ^ 

24m 
+0( 7 2 



(1 + 3C 2 ) + ^(1-C 2 )7 



4m 



The strongly attractive limit from our parametrization 
formula is 



£Gs{n, C,7 -> -°°) 



(1-C)7 2 



8m 

h 2 ir 2 n 2 1 



which has little difference comparing with Batchelor's 
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derivation 24 1 



x(l + 15C 3 



C)7 2 - 
-3C 



3C 



96m 

* ■ 0(1). 



For emphasizing these limits, we show egs{ti, C> 1)/Ef 
as functions of £ and 7 in Fig. [2] At £ = 0, 
£GS (i) 0, j)/Ep behaves like —7 at weak attractive limit 
and — 7 2 at strong attractive one, while at £ = 1, 
£gs(?i, 1, j)/Ep goes as a constant. The expression at 
strong interaction limit differs from the corresponding 
repulsive case due to the attractive nature of the inter- 
action. Furthermore eQs(n,£,~f)/Ep behaves like £ 2 a t 
weak attractive limit similar to the corresponding repul- 
sive ID expression, however it behaves like £ at strong 
attractive limit. 

Furthermore, we would like to examine some other 
physical quantities calculated from Eq. ([7|) in comparison 
to those calculated within the exact Bethe-Ansatz equa- 
tions to obtain more credibility for our parametrization 
formula. For this purpose, we take the derivative of the 
ground state energy with respect to the magnetization 
to calculate the magnetic field in equilibrium condition 
which is h(n,s,j) = d[neGs( n , C l)]/ds [25|. The mag- 
netization vanishes when the field becomes smaller than 
the critical value h c , a term associated with the spin en- 
ergy gap in the attractive case (Note that it vanishes in 
the repulsive case.). We next display the exact magne- 
tization as a function of the magnetic field solved nu- 
merically from the coupled integro-differential equations 
based on Eqs. ©-© in comparison to those derived from 
the parametrization formula. 

For given 7 = —10, the exact numerical calculation 
gives h c — 6.097 for n — 0.5, while from the parametriza- 
tion formula we get h c = 6.049 for the same density as 
shown in Fig. [3] Moreover, h c = 24.388 for n = 1.0, by 
calculating exact formula and we have h c — 24.182 from 
parametrization formula. Surprisingly, a finite magneti- 
zation appears only if the magnetic field larger than h c . 
The magnetization saturates at s — n/2 when h > h s . 
More precisely, for 7 = —10, the exact solution of h s is 
8.410 for n = 0.5, and 33.643 for n — 1.0, which are 
very close to values obtained from the parametrization 
formula that we get h s — 8.440 and 33.759, respectively. 

After comparing the exact solutions to the 
parametrization formula, we convince that the 
parametrization formula gives satisfying simulations on 



the ground state energy and the derivative quantities 
like the chemical potential and the magnetic field in a 
wide range of parameters. 



III. SPIN-DENSITY-FUNCTIONAL THEORY 
AND THOMAS-FERMI APPROXIMATION 



The ground-state spin densities n a {z) can be calcu- 
lated by spin-density-functional theory (SDFT) solving 
self-consistently the Kohn-Sham equations 



Ira dz 2 



Va,<y{z) = e a ,cripa,a(z) (12) 



with effective potential (2) = ^h°^ bvK*) + 

Vxc [ n a-]( z ) + rauj 2 z 2 /2. The ground state density is de- 
termined by the closure 



■cor 



(13) 



a=l 



Here the sum runs over the occupied orbitals and the de- 



la) 



generacy factors r„ satisfy the sum rule J2 a 
The spin-dependent effective Kohn-Sham (KS) potential 
as usual is composed by the mean-field term V^j = 
g\T>ns{z), the exchange-correlation potential defined as 
the functional derivative of the exchange-correlation en- 
ergy E xc [n a ] evaluated at the ground-state density pro- 
file, Vxc* 1 = SE- xic [n CT ]/dn a -(z)\ GS , and the external poten- 
tial, respectively. 

In order to calculate n a (z), E xc needs to be approxi- 
mated. Here we take the local-spin-density approxima- 
tion (LSDA) which is known to provide an excellent de- 
scription of the ground-state properties of a variety of in- 
homogeneous systems [3] to calculate. In the following 
we employ an LSDA functional based on the parametriza- 
tion results of Eq. |(7J) for the exchange-correlation poten- 
tial, 



£ X cK] = J dzn(z)e h x r(nXn)\ n ^ n{z) , ( _> c{z y (14) 

The exchange-correlation energy e^° m of the homoge- 
neous Gaudin-Yang model is defined as 

^° m ( n , C, 7) = £Gs(n, C, 7) ~ n(n, C) - En(n, C, 7) ■ 

Here £h(«,Cj7) = T^~7(l — C 2 ) comes from the contri- 
bution of the Hartree-Fock mean field. 
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Following the above points, we write down conve- 
niently the Kohn-Sham potential as, 

1 



muJuZ + /j,[n, s](z) — no[n, s](z) 



where fi[n,s](z) = 
is the chemical 
d[ns GS (n,(,j)}/ds 



±^{h[n,s](z)-h [n,s](z)} , (15) 

= d[neGs(n,C,i)]/9n\n-m(z),s-^a{z) 
potential and h[n,s](z) = 



n^n{z), s ^s{z) the magnetic field 
of the homogeneous system evaluated at the local 
potentials, respectively. [io[n,s](z) and h [n,s](z) are 
the corresponding chemical potential and magnetic field 
of the noninteracting system, which are given by 



fj, [n,s](z) = 
h [n,s](z) = 



8m 
H 2 ir 2 



[n 2 (z)+As 2 (z)] 
n(z)s(z). 



Here, n(z) = n-\{z) + ni(z) is the total density distri- 
bution, and s(z) = \n^{z) — ni(z)]/2 is the spin density 
distribution or the local magnetization. 

From the discussions in Sec. (|IT]) . we know that our 
spin-density-functional theory with the reference system 
based on Bethe-ansatz equations has correctly incorpo- 
rated the Luther-Emery-liquid nature and the spin en- 
ergy gap of the homogeneous system 27, 0|. The in- 
homogeneity nature induced by the harmonic trap will 
be incorporated into the exchange-correlation potential 
through the local spin-density approximation. Hereafter 
we refer to this method as BALSDA. 

Another simplified method to describe this inhomo- 
geneous system is to resort to a local density approxi- 
mation (LDA) also for the noninteracting kinetic energy 
function, which is a Thomas-Fermi approximation (TFA) 
similar to the TFA but implementing the exchange- 
correlation energy. We refer to this approach as the TFA 
in the following. Within the LDA one assumes that the 
chemical potential of the trapped system is given by the 
sum of the spin-dependent local chemical potential taken 
to be the chemical potential of the uniform system at the 
corresponding density, and the external potential. The 
ground-state density profile is obtained in the spirit of 
the local equilibrium condition, 



(16) 



derived by directly minimizing the total energy functional 
of the system. The constants /i° are fixed by the normal- 
ization condition N a = J dzn a (z), and /j,^ om [n^, njj(z), 



the corresponding density dependent chemical poten- 
tials of the homogeneous system is derived from the 
parametrization formula in Eq. ((7]) evaluated at the local 
densities of n-\{z) and ni(z). Within the LDA, a dimen- 
sionless characteristic parameter r/ — Na 2 D /a 2 can be 
defined in the harmonic trap. The systems hold the same 
density profile if the characteristic parameter is the same 
irrespective of their different number of particles and dif- 
ferent harmonic oscillators [lp| . The parameter can be 
also expressed as, r\ — AN/ A 2 with A = gi£,/(a\\hid\\) if 
we keep in mind that 7 = —2/ (naio)- T) ^> 1 corresponds 
to weak coupling while r\ <C 1 corresponds to the strong 
interacting regime. 

IV. NUMERICAL RESULTS 

In this section we present the numerical results ob- 
tained from the self-consistent solution of Eqs. (|12 |l - (fl"5|) 
using our accurate parametrization formula from Pade 
approximant. We test the TFA approach and compare 
its results with our full BALSDA to justify how advance 
physics we have within BALSDA. We proceed to illus- 
trate our main numerical results, which are summarized 
in Figs. 4-5. 

First of all, we show the ground state site occupation 
for a Fermi gas trapped in a harmonic potential at differ- 
ent polarizations in Fig. 2J Here A = — 2 which implies 
rj = A, corresponding to weak coupling regime. For the 
harmonically confined system, we define the average po- 
larization as P = (A| — A|)/(Af + Ni), which is changed 
keeping always a constant number of spin-up atoms JVj 
and decreasing N±. P = means an unpolarized system 
and P — 100% corresponds to a trivial fully polarized 
case of Af noninteracting fermions. The density profiles 
of spin-up and spin-down atoms show Aj and max- 
ima, respectively 11], which reflects in the total density 
profile with N± maxima and in the spin one with N± 
minimum. 

In Fig. IDJa), the total charge density profile is shown 
as a function of z. The 'hat '-structure is formed in 
the bulk region by paired minority species of spin-down 
and spin-up atoms, manifested by N± peaks inside and 
partially polarized fermions. Higher polarization gives 
smaller 'hat'-structure. By increasing the polarization P 
but keeping iVj constant, the range of the charge density 
becomes a little larger because the attractive interaction 
becomes weaker due to the decreasing spin down atoms. 
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(b) 
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FIG. 4: (color online) The charge density (a) and spin density 
(b) profiles as a function of z for different polarization P at 
A = —2. The polarization is increased from 3% to 90% by 
keeping always a constant number of spin-up atoms, N] = 20 
and decreasing N\, from 19 to 1. 



Orso 14[ and Liu et al. la, [29( , after analyzing the phase 
diagram and calculating the order parameter, concluded 
that the state at the center of trap gives a two-shell struc- 
ture of the FFLO-type, which has a spatially oscillatory 
behavior in the superconducting correlation function and 
exists in any nonzero polarization and attractive inter- 
action. In Fig. 0Jb), we show the spin density profile 
defined as the difference between spin-up and spin-down 
density, namely the local magnetization as a function of 
z. By increasing the P value, an oscillating structure of 
the spin density wave appears in the bulk of the system, 
due to the attractive interaction between the majority 
and minority species. As a result, a quasi-flat region ap- 
pears in the bulk part of the spin density profile with 
coexisted spin-up and spin-down atoms, surrounded by 
fully polarized spin-up atoms only. The core of coexist- 
ing spin-up and spin-down atoms becomes smaller and 
smaller for higher value of P. We would like to point out 
that the local magnetization distribution as a function 
of position is directl y m easurable through phase-contrast 
imaging techniques [30(. 



In Fig. [5l we show the spin-up, spin-down density pro- 
files, and furthermore, the local magnetization distribu- 
tion at A = —12 with fixed total atom number of N = 36 
but varying spin-up and spin-down atoms. This implies 
rj = 1, corresponding to an intermediate coupling regime. 
From Figs. (Ha), (b) and (c), it is clear that the ampli- 
tude of the oscillations increases with decreasing P. We 
present results obtained by the TFA approach together 
with BALSDA to justify quantum effects. The TFA re- 
sults are shown by solid lines between the oscillating den- 
sity profiles. For P = 0.56, see Fig.[5£a), the TFA results 
give a good representation of the actual density profile, 
except at the edges of the cloud. We stress here, without 
including the exchange-correlation effect, the TFA gives 
bad overall shape for the density profiles. 

In Fig. EJa), a high polarization case of P = 0.56 
is presented for the spin-up density, spin-down density 
and the spin-density profiles. It is clear that the curve 
envelop produces an oscillating structure particularly in 
Figs-E^b) and[5^c). A two-shell structure is formed with 
a partially polarized FFLO superfluid in the center of the 
trap Q and a fully polarized normal state of excess spin- 
up atoms at the edges. In Fig. EJb), a polarization case 
of P = 0.11 is illustrated with BALSDA together with 
TFA, where a different two-shell structure is formed with 
a partially polarized superfluid at the trap center and a 
fully paired BCS superfluid state at the edges, opposite 
to the 3D counterpart of the BCS paired phase in the 
core. In the present system, we observe that there is 
one critical point which is P c = 0.12 between the two 
above-mentioned different phase separations which coin- 
cides with the calculations given by Refs. BQ. We 
want to stress here, from Fig.[5]Ja) to Fig.^b), due to the 
stronger attractive interaction when the size of the cloud 
shrinks considerably by decreasing P, that the oscillation 
structure due to correlation effects become more impor- 
tant and the LDA result becomes worse. In Fig. 02c), an 
even smaller polarization case of P = 0.06 is illustrated 
with BALSDA against TFA results. The fully paired su- 
perfluid state becomes larger in the wings with such a 
smaller polarization. The different phases in the two- 
shell structure of standard BCS state outside of the cen- 
ter and the FFLO-type state in the bulk, are proofed to 
be a smooth second-order transformation [29!| . In all the 
above calculations with weak and intermediate coupling 
regimes, we find the TFA results give a good overall shape 
to the microscopic calculations based on the BALSDA. 
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FIG. 5: (color online) The spin-up, spin-down and spin den- 
sity profiles of fixed atom number N = 36 for r) = 1. We show 
two different results based on the BALSDA (dotted lines) and 
the TFA (solid line between the oscillatory BALSDA results) . 
(a) For P — 0.56, a two-shell structure is formed with a par- 
tially polarized FFLO-type superfluid of distinctive oscilla- 
tory character in the center of the trap and a fully polarized 
state at the edges, (b) For P — 0.11, a two-shell structure is 
formed with a partially polarized superfluid in the center of 
the trap and a fully paired superfluid state at the edges, (c) 
For P = 0.06, a smaller polarization gives larger BCS regime 
in the wings. In all the calculations, we show the TFA results 
give a good overall shape between the oscillating BALSDA 
ones. 

But importantly, the TFA could not produce oscillation 
structure due to correlation effects in the quantum sys- 
tem. 



V. CONCLUSIONS 

In conclusion, we have presented two alternative 
methods to study the imbalanced two-component Fermi 
atomic gases of attractive interactions in Q1D harmonic 
traps. One is the spin-density-functional theory based on 
Bethe-ansatz solutions, which properly incorporates the 
Luther-Emery nature and the spin energy gap of the ho- 
mogeneous system. Another is the TFA which takes a lo- 
cal density approximation for the noninteracting kinetic 
energy but including the exchange-correlation energy. 
We apply these two methods to the system of finite size 
trapped by the harmonic confinement. At weak coupling 
strength of any polarization or intermediate coupling of 
large polarization, a two-shell structure is obtained with 
a partially polarized pairs of FFLO-type state in the core 
and a fully polarized fermions in the wings. At interme- 
diate coupling regime of small polarization, a two-shell 
structure is formed of partially polarized pairs of FFLO- 
type state in the bulk and a fully paired BCS state at the 
edges which is different from the 3D case. In the current 
experimental set-up, imbalanced interacting Fermi gases 
can be induced by a radio-frequency sweep in the opti- 
cal lattice system, and in the on-going experiments, the 
two-shell structure described here should be observable 
in the system of thousands of Q1D tubes formed in the 
quasi-two-dimensional optical lattices with each of the 
tube containing up to 100 atoms measured by absorption 
imaging techniques, while the direct test of FFLO states 
remains as a challenge to experimental community [311 ]. 

It would be of interest to develop the present scheme to 
study dynamical phenomena in these strongly correlated 
gases using time-dependent DFT and/or current-DFT, 
instead of resorting to the inhomogencous Tomonaga- 
Luttinger liquid model. From a more formal DFT view- 
point, a functional better than the one presented in 
Eq. {T3J) is desirable and necessary to deal with the strong 
coupling regime. 
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